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ABSTRACT 

We present the results of an attempt to adapt the distribution function formalism to 
characterize large-scale structures like clusters of galaxies that form in a cosmological N- 
body simulation. While galaxy clusters are systems that are not strictly in equilibrium, 
we show that their evolution can nevertheless be studied using a physically motivated 
extension of the language of equilibrium stellar dynamics. Restricting our analysis to 
the virialized region, a prescription to limit the accessible phase-space is presented, 
which permits the construction of both the isotropic and the anisotropic distribution 
functions /(£) and f(£,L). The method is applied to models extracted from a cata- 
logue of simulated clusters. Clusters evolved in open and flat background cosmologies 
are followed during the course of their evolution, and are found to transit through a 
sequence of what we define as 'quasi-equilibrium' states. An interesting feature is that 
the computed f(£) is well fit by an exponential form. We conclude that the dynamical 
evolution of a cluster, undergoing relaxation punctuated by interactions and violent 
mergers with consequent energy-exchange, can be studied both in a qualitative and 
quantitative fashion by following the time evolution of /(£)■ 
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1 INTRODUCTION AND MOTIVATION 

In this paper we present the results of an attempt to con- 
struct distribution functions for clusters of galaxies that 
form in N-body simulations. To permit the description of 
systems like clusters that are continually evolving at the 
present epoch we propose a possible extension of the distri- 
bution function formalism of equilibrium stellar dynamics. 
The prescription developed is applied to simulated clusters 
for which the isotropic and anisotropic distribution func- 
tions are computed. We then use the time evolution of the 
distribvj+i Hit -fa iilLiuu Lu ycL an iiisiuhL uiLu the details uf Lin 
relaxat 



dynamical state of clusters formed in N-body simulations, 
(ii) studies of the formation and evolution of elliptical galax- 
ies in an attempt to understand the physical origin of the 
universal R 1 ^ law and (iii) the evolution of stellar systems 
studied in phase space. 

Given the standard paradigm of Cold Dark Matter 
(CDM) dominated hierarchical structure formation scenar- 
ios, the formation, internal structure and evolution of dark 
halos and density profiles has been studied in detail for 
a range of initial input power spec tra evolved in a slew 



of background cosmological models (Goldreich 1984 



pmmm and ewg 

Th'e motivation for this work is twofold: (i) to develop a 
mathematical formalism to characterize the dynamical evo- 
lution of 'stellar' systems that are marginally out of equi- 
librium, which can be applied to numerical cluster models 
and (ii) to gain insight into the dynamical evolution and 
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man & Shaham 1985; West, Dekel, & Oemler 1987; Crone 



Hoff- 



Evrard, fc Richstone 1994|; Richstone, Loeb, fc Turner 199S 



prop erthes ot clusters m general witii as special emphasis 
on the origin of density profiled and the Use of Clusters as 
cosmological probes. 

Previous studies that are relevant to our work have been 
done in three somewhat distinct contexts: (i) studies of the 



Cen, Gnedin, fc Ostriker 1993t |Klypin et al. 1993fc pen 199j ; 
Navarro fc White 1994} podds 1995[ ). A 'universal density 
profile' has been found to be a good fit over a two-decade 
range of scales for dark halos th at form in N-body realiza- 
tions of the cosmological model (Navarro, Frenk, & White 



199C White 1996 ). These authors and several other gr oups 



"( Pole fc Lacey 1996| |Tormen, Bouchet, fc White 1996|) find 
that for a range of cluster-halo shapes, sizes and other bulk 
properties, the profile depends primarily on the epoch of 
formation and only rather weakly on the details of the cos- 
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mological context and the initial power s pectrum of flue- 
White (1996| ) find that their 



tuations. 



Navarro, Frenk, 



spherically averaged density profiles can be fit by scaling a 
'universal' profile (henceforth referred to as the NFW pro- 
file) which goes as p oc r _1 at small radii and steepens to 



p oc r at large radii. Subsequently, Cole & Lacey (199C ) 



on examining the internal structure of dark matter halos 
formed in scale-free hierarchical models for a range of initial 
power spectra, also find that the density profile is well-fi t 
by both the NFW profile and the analytic Hernquist (199C ) 



profile which falls off as p oc r~ at large radii. Part of 
the aim of the present investigation is to explore the inter- 
play between the initial conditions and subsequent relax- 
ation processes in clusters and to understand the possible 
physical origin of these density profiles. 

This approach is somewhat analogous to the det ailed 
studies of the origin of the R}l 4 law in elliptical ga l axies (|van 
Albadal 19821 ICarlbcrg, Lake , fc Norman 198(1 |Londrillo 



Messin;,, & Stiavelli 1991; Hjorth & Madsen 1991). Moving 



over to the analogy with the formation of elliptical galaxies, 
we briefly review some scenarios that have been put forward 
to understand the observed universal R 1 ^ 4 profile. The phys- 
ical system in this case essentially consists of stellar orbits in 
a potential which deepens during the course of evolution via 
violent relaxation that occurs as a consequence of strong 
collective potential fluctuations ( Lynden-Bell 1967 ). Dur- 
ing the course of evolution following the violent phase, any 
remaining positive energy particles slowly diffuse out and 
eventually leave the system. In the case where the resultant 
stellar system has a deep potential well the R 1 ^ 4 surface- 
brightn ess pro file arises as a natural consequence ( Hjorth & 
Madser 1991 ). Such a final configuration can be produced 
from a cold dissipationless collapse scenario leading to stel- 
lar orbits that are preferentially radially anisotropic as well 
as from warmer initial conditions with dissipation. 

While universality intuitively suggests wiping out of the 
memory of initial conditions and similarities in the post - 
collapse evolution, Londrillo, Messina, & Stiavelli (1991) 



conclude from their study of dissipationless galaxy forma- 
tion that the phase-space properties of the final equilibrium 
state do in fact depend strongly on the details of the initial 
conditions. They find that the relaxation to an R 1 ' 4 profile 
is characterized by two phases: (i) an early phase dominated 
by strong potential and density fluctuations in the central 
region which results in a core-halo structure and a long tail 
in both the density and energy distribution and (ii) a longer- 
lived later phase wherein the system evolves mainly due to 
the effect of small-scale diffusion which smoothes the con- 
nection between the core and the halo populations. Some 
memory of the early phase is retained, since the fluctuat- 
ing forces and the compact core necessary to diffuse the low 
angular momentum orbits in the later phase are necessarily 
remnants of the initial collapse dynamics; i.e. the region of 
phase space that is available for diffusion into for the system 
during the late phase is delinea ted during th e early phase. 
Also in dissipational scenarios (Barnes 1996) which are ex- 
pected to introduce irreversibility, especially during mergers, 
some memory of the initial conditions seems to be retained. 

Finally, we place this work in the context of studies of 
the evolution of stellar systems in phase space. The structure 
of phase s pace and its a ccessibility for R 1 ^ 4 galaxies was ex- 
plored by Binney (1982 ) who found that the distribution of 



the number of stars of a given energy is approximated rather 
well by an exponential formula. It was argued that if ellip- 
tical galaxies formed via dissipationless processes, the expo- 
nential form found at late times for the distribution function 
was probably set up during the epoch of formation, yielding 
a scale-free form for the density profile. Therefore, the phase 
space corresponding to the observed surface brightness pro- 
files of elliptical galaxies are consistent with having been 
demarcated during the initial collapse of a density fluctua- 
tion t hat subsequently gives rise to the g alaxy. In an another 
study, Hernquist, Spergel, fc Heyl (1993 ) attempted to quan- 



tify the importance of mergers in the origin and structure of 
early-type galaxies by exploring the phase-space properties 
of merger remnants and comparing them to simple models 
of elliptical galaxies. 



In a more recent analysis, Vbglis (1994 ) presents an ana- 
lytic model fit to the distribution function of a nearly spher- 
ical, anisotropic galaxy-scale system. Starting with cosmo- 
logically consistent initial conditions, he finds a dichotomy 
of states, a core population dominated by low angular- 
momentum tightly bound particles and a less bound halo 
population. The energy distribution of the core population 
is nearly isothermal whereas the halo population consists 
of particles on radial orbits with correlated energy and an- 
gular momentum. Following the time evolution of a radius 
containing a specific fraction of the mass, he finds that the 
radial oscillations decay quickly after collapse preceding the 
onset of relaxation during which a new equilibrium is estab- 
lished. 

Despite some essential physical differences between stel- 
lar systems and clusters, we synthesize the approaches de- 
scribed above and demonstrate that they can be adapted 
and successfully applied to the case of galaxy clusters that 
form in realizations of cosmological N-body simulations. 
Specifically, we develop in the present paper a framework 
that can be applied to characterizing what we shall define 
as quasi- equilibrium evolution. In Section 2 we outline our 
conceptual model of a cluster, review the equilibrium distri- 
bution function formalism and present our prescription to 
extend it to quasi-equilibrium. The details of the N-body 
simulations from which the cluster catalogs that we study 
were extracted are described in Section 3. In Section 4 we 
apply the method to the construction of f(£) and f(£,L) 
for slices from the simulations. We present our results and 
summarize what can be learnt from a first set of diagnos- 
tics developed to describe the quasi-equilibrium evolutionary 
states of the cluster in Section 5. 



2 THE DISTRIBUTION FUNCTION 
FORMALISM 

2.1 Characterizing a cluster 

A cluster is a composite system with diverse components: 
collisionless dark-matter particles, collisional baryonic gas 
and bound stellar systems (the galaxies). The rudimentary 
physical picture of the cluster that we find useful, is one of 
a system with two natural length scales delineated by the 
dynamics: the virial radius r v i r , defined to be the radius en- 
closing an overdensity of « 180 times the critical density, 
which is motivated by the collapse of the simple spherical 
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top-hat model (Peebles 198C) and the turnaround radius 



rta, defined to be the radius at which the cluster separates 
from the general cosmological expansion (envelope with zero 
relative velocity). There is on-going infall inward of r ta , and 
therefore the system can never really be treated as one evolv- 
ing in isolation. 

It is instructive to clarify at this juncture some of the 
terms and their definitions that will be employed in our 
analysis to describe the dynamical state of the system. A 
virialized system is one for which the virial theorem holds 
(d 2 Iij/dt 2 = 0; where Iij is the moment of inertia tensor) 
and therefore there exist no systematic motions, neither ex- 
pansion nor contraction; a stationary system is one for which 
there is no time variation of the distribution function /, i.e. 
df jdt — 0, and a relaxed system is one for which memory 
of the initial conditions in phase-space has been erased. For 
a cluster-scale collisionless system, two-body relaxation pro- 
cesses are unimportant since the time-scale relevant for two- 
body processes far exceeds the dynamical time, although in 
the case of a simulated cluster, there arise some purely nu- 
merical two-body effects which we discuss in more detail in 
Section 3. While collective energy exchange processes might 
occur during the course of dynamical evolution, they might 



where df/dt is the total derivative along the exact trajec- 
tory of the particle moving in the force field V5> and C is 
the collisional term, which defines the changes in / due to 
collisions between particles. In conjunction with the Poisson 
equation, 



V 2 $ = A-nGp, 

and the definition of the density profile, 
p(r,t) = / ,f(r,v,t)d 3 v, 



(2) 



(3) 



equations (1), (2) and (3) constitute a complete set of 
self-consistent evolution equations that describe a self- 
gravitating system. For a stellar system that is strictly col- 
lisionless and gravitating, C = in equation (1), the solu- 
tions to the above are the stationary states, and the parti- 
cle trajectories are the 'characteristics'. Additionally, while 
equation (2) is valid only for a self-gravitating system, the 
potential defined in equation (1) could also include the con- 
tribution from an external potential. 

We review below some of the basic definitions that will 
be used throughout the paper; a more c omprehensive treat- 



not alw ays lead to a relaxed final state and also a virialized 



ment of the formalism can be found in Binney & Tremainc 
(19871) . The notational conventions followed here are also as 



system need not be relaxed. This is so because virialization 
is defined in terms of the kinetic and potential energies of a 
specific configuration of the system at a given epoch, while 
relaxation requires physical processes that affect the evolu- 
tion to be taken into account and is hence a statement about 
the integrated dynamical history. 

In this study we focus on the collisionless component. 
We define a collisionless equilibrium system to be one that 
is virialized, stationary and relaxed with no on-going energy 
exchange. Clearly, this is not the case for a cluster of galax- 
ies, wherein departures from equilibrium necessarily occur 
due to the existence of an infall region and departure from 
stationarity which is a consequence of fluctuations induced 
in the potential in response to mergers and secondary infall. 
Nevertheless, there does seem to exist a virialized central 
region and observed clusters seem not to be far from equi- 
librium as inferred from their 'measured' density profiles and 
velocity dispersions. Therefore, the need arises to develop a 
phase-space description that can incorporate small depar- 
tures from equilibrium, which we attempt below. To first 
approximation, we model a cluster in this treatment as a 
spherically symmetric and non-rotating system. 



2.2 The equilibrium distribution function 
formalism 

For a statistical description of a system with a large num- 
ber of particles, it is convenient to define the distribution 
function (DF) (the phase-space mass density) f(r,v,t). 
A given configuration of the system is then specified by 
f(r, v, t)d 3 xd' i v - the mass of particles having positions 
in the infinitesimal volume element d 3 r, with velocities in 
the range d 3 v in the (r, v) phase space. The DF satisfies the 
Boltzmann equation, 



in Binney fc Tremaine (1987 ) except for the introduction of 
the function h. 

The binding energy per unit mass £ is defined as, 

£ = *(r) - y, (4) 

and the angular momentum per unit mass L is given by, 

L = |r x v|. (5) 

The DF for a spherically symmetric, non-rotating stellar sys- 
tem can only depend on these two quantities. Dependence 
on both £ and L gives rise to an anisotropic velocity dis- 
tribution whereas dependence on £ only yields an isotropic 
distribution of velocities. The differential energy distribu- 
tion {dM(£)/d£)d£ is the total mass of the system with 
binding energy within [£ , £ + d£] which is then explicitly 
defined to be, 



K £)d£ = ^§±d£ 



f{£)d\d\, 



(6) 

integrated over the volume element AV(£) in phase-space 
with energy in the specified range [£, £ + d£ ]. Since we are 
dealing with spherically symmetric systems we can write this 
integral in terms of r and v. Furthermore, we can take f(£) 
out of the integral and use eq. (4) to arrive at: 
rr m (E) 

h(£)d£ = I6iv 2 f{£)d£ / r 2 dry/ 2 - £). (7) 

Jo 

This can be written as 

h(£)dS = f(S)g(S)dS, (8) 

where g(£ ) the density of available states in phase space of 
energy £, 

g(£) = 16tt 2 / y/2 (tf (r) - B) r 2 dr. (9) 



df df _ _ 8/ 



(i) 



This integral is strictly defined only for 9(r) > £ . Here 
r m (£) is defined to be the radius where ty(r m ) = £. Thus, 
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we obtain an operational definition for /(£), 



f{£) 



(10) 



Given a potential *l/(r), both h(£ ) and g(£ ) are well-defined 
for stellar systems in equilibrium. We use the definition in 
equation (9) to construct the DF in our analysis. The DF 
could alternatively be derived and hence constructed by ei- 
ther explicitly counting available and energetically accessible 
cells in phase space or from the density profile of the sys- 



tem pjr) using an inversion method devised by Eddington 



(1916 ^jQuoting th e result of Eddington's inversion formula 
from Binney (198^ ), 



d 2 p = f r 2 

d* 2 ycM 



d* 



+ 



1 dp 



d* 2 y/W^£ ^d* 



dr 2 



4-Kpr' 
M{r) 



(11) 



(12) 



where M(r) is the total mass enclosed within radius r of the 
system. 

For a spherically symmetric system with an anisotropic 
velocity distribution, the DF f(£ , L) is analogously defined 
to be, 



/(£, L) 



h(£, L) 



(13) 



The density of states in phase-space g(£ , L) is now a hyper- 
surface of energy £ and angular momentum L. Furthermore 



g(£, L) = 8ty 2 LT t (£, L), 

where T r is the radial period defined as, 

T t {£,L) 1 '''* 



L- 



(14) 



(15) 



the range of integration being from the pericenter of the or- 
bit, n, to the apocenter, r^- The radial period T T (£, L) is 
generally evaluated numerically, but an approximately ana- 
lytic expression can be obtained for a few special choices of 
the potential; one of them being the isochrone potential, for 
which the following simple expression holds, 

m, v> = UK- m 

Furthermore, T r (£,L) has been demonstrated to be almost 
independent of L for diffe rent potentia ls (is ochrone and 
Keplerian in particular) by Voglis (1994) and Lynden-Bell 
(19960 



2.3 



Adapting the DF formalism to 
quasi-equilibrium 



domain. Indeed, most studies of the dynamics of clusters of 
galaxies implicitly assume that an equilibrium description is 
an adequate starting point. 

Below we first describe the inherent limitations that 
such an approach is necessarily subject to, but we then use 
precisely these restrictions to formulate a physically moti- 
vated and self-consistent extension of the equilibrium for- 
malism. The key to the extension arises from the recognition 
of the important distinguishing features in the case of clus- 
ters: the finite time available for evolution and the spatial 
extent of a cluster. 

Examining the relevant time-scales: firstly, the typical 
time that it takes for a cluster to be virialized out to the 
Abell radius is of the order of the Hubble time which in it- 
self precludes the cluster from reaching an equilibrium con- 
figuration. Secondly, the disturbances from continuous infall 
and the more violent interactions constantly drive the sys- 
tem away from equilibrium whereas the collective energy 
exchange processes tend to drive the system toward equilib- 
rium. The fact that the characteristic time-scale on which 
the collective energy exchange process operate is relatively 
short means that the system is capable of quickly recouping 
from or adapting to a disturbance. Thus, it is possible that 
the system might come close to equilibrium fairly quickly 
after such events. In other words, it may be possible to de- 
couple the system into two components: a fast component 
and a slowly evolving component, akin to the adiabatic ap- 
proximation. We may therefore view the system as evolving 
from one quasi- equilibrium state to another on a time scale 
that is roughly the 'disturbance time scale'. 

Having thus argued that the system might frequently 
be close to a quasi-equilibrium state, we need to delineate 
more precisely the specific region to which such a treatment 
is most suitably applicable. Since the use of the equilibrium 
formalism necessarily requires the system to be virialized, a 
natural choice would be to study a typical virialized region; 
clearly, the infall region cannot be studied via this approach. 
Also, once a radially infalling particle has entered the virial- 
ized region its kinetic energy gets distributed to other par- 
ticles and it is more likely to get bound inside the virialized 
region rather than leave the system. Moreover, even the few 
particles that are capable of leaving the region are more than 
balanced in the net by the infalling ones. 

In summary, the finite time and extent arguments cou- 
pled with the expected short relaxation time scale encourage 
us to develop the notion of quasi-equilibrium states further. 
However, the eventual usefulness and indeed validity of the 
approach should be judged not on these arguments alone 
but on the results presented in Sections 4 and 5. We reiter- 
ate that our proposed adapted equilibrium formalism is only 
applicable to a restricted region in phase space. 

There have been several previous studies that have at- 
tempted to limit the available phase space for elliptical 



The eq lilibrium formalism outlined above does not strictly 



galaxies ( Stiavclli & Bertin 1987; Fremaine 1987 ; Merritt 



apply t ) clusters of galaxies as these systems are constantly 



evolving. Additionally, the equilibrium formalism as cur- 
rently posed cannot accommodate the existence of positive 
energy particles that are sometimes observed in numerical 
cluster simulations (see plots in Section 5). On the other 
hand, it appears to be sensible to attempt using the equi- 
librium formalism as a guide to probe the non-equilibrium 



Tremaine, fc Johnstone 1989). For instance, Spergel & Hern 
quist (1992j ) suggested a prescription to limit the macro- 



scopic states that are accessible to individual particles in a 
collisionless system approaching equilibrium by approximat- 
ing violent dynamical processes en-route to relaxation as a 
sequence of discrete scattering events. They also introduced 
a monotonically increasing non-equilibrium entropy by im- 
posing additional constraints on the available phase-space. 
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In our approach we make no a priori assumptions about 
the nature of energy-exchange processes and do not attempt 
to define any generalized entropy for the system. Rather we 
shall argue that it is the usual definitions (8) and (13) of g(£ ) 
and g(£,L) that need to be modified in a consistent fash- 
ion to enable the description of systems evolving marginally 
away from equilibrium. 



2.3.1 Computing g(£) for non- equilibrium systems 

As described in Section 2.2, the phase-space density g(£ ) 
is valid only for negative energies (positive £ , i.e. strictly 
bound particles in an equilibrium system), and diverges for 
£ — > due to the blowing up of the volume element in 
physical space. For instance, consider a potential ^(r) with 
a radial dependence such that, 



*(r) = r~ 



< a < 1. 



(17) 



In this case r m (£) ~ £~ a and so g(£ ) diverges as £ v for 
r\ — a 2 /2—3a, i.e., —5/2 < r\ < 0. Requiring f(£) to be finite 
at the escape energy and the mass to be finite simultaneously 
implies that 77 > — 1, in turn restricting < a < 3 — Vo 
which for the self-gravitating case (a = 1) means that the 
asymptotic behavior of g{£) cannot be compensated by a 
similar divergence of h(£). Furthermore, g(£) is infinite for 
all £ < 0, as the accessible phase-space is infinite for these 
energies. As discussed above, we propose to limit the acces- 
sible phase-space due to the finite time available. We do this 
in practice by truncating g(£ ) at a fixed scale r max . As we 
show below, this will remove the divergence and render g(£) 
well-behaved for all £ . 

Physically, such a truncation amounts to the assump- 
tion that for most of the bound particles that constitute the 
system, the finite time elapsed since the initial collapse nec- 
essarily implies that the probability of occupation of states 
falls off steeply beyond a fiducial volume in phase space 
which is dictated by the threshold distance that a parti- 
cle with typical velocities can traverse in the available time 
interval. Energy-exchange and hence marginal departures 
from equilibrium of the cluster arise both from the existence 
of the extended infall region as well as due to the occurence 
of frequent mergers and interactions as an integral part of 
the relaxation process. Since, in this analysis, we are primar- 
ily interested in probing and characterizing the evolution of 
the central regions and the relaxation effects therein, we 
choose to truncate the system at r max = riso = r v ir, the 
radius within which the system is expected to be virialized, 
and is indeed found to be for several different power spectra 
(Cole & Lacey 1996). This prescription is independent of the 



functional form for the potential and can in general be ap- 
plied to any physically meaningful potential for a spherically 
symmetric system. 

Specifically, consider a stellar system consisting of par- 
ticles confined in a potential ^(r), modelled explicitly as the 
Hernquist potential, 



*(r) 



#0 



1 + 



(18) 



In the equilibrium context, the limits of integration for g(£) 
range from to r m (£) defined to be the radius at which 
*(r«) = £■ 



-5x10 s 





[km 2 s~ 2 



Figure 1. The prescription for truncation: we illustrate the pre- 
scription we have used to define r max , the truncation radius, for 
the Hernquist potential. The solid curve is our choice for r max , 
which now stands well-defined for all C. 



From the plot in Fig. 1, we see clearly that the phase- 
space volume element is ill-defined as £ — > + . By setting 
r max as the maximal allowed value for r m {£) in eq. (8), i.e., 

1 £ 



,(£) = Min 



-1) 



(19) 



we truncate the system in physical space, which automat- 
ically translates into a restriction of available phase-space. 
In this way we recover a well-behaved g(£ ) that is valid for 
all £. 

Thus, using only the assumption that the number of 
particles within the boundary r max is roughly constant, we 
modify the definition of g(E) to allow for quasi-equilibrium 
states (including the temporary existence of positive energy 
particles). With this modification, most of the equilibrium 
formalism then automatically carries over for spherically 
symmetric systems. 



2.3.2 Computing g(£, L) for non- equilibrium systems 

The construction of g(£,L) is analogous to the above, but 
we now need to restrict a two-dimensional hypersurface. The 
key observation here is that in addition to the restriction in 
?"m ax , the orbital structure needs to be constrained corre- 
spondingly in a consistent fashion. We proceed as follows, 
by first assuming that g(£,L) is separable, 



g(S,L) = Lg(£), 



(20) 



motivated by the functional form of the expressio n for 
g(£,L) (eq. 13), and the observation of Voglis (1994) that 
the radial period T r (£,L) is only a weak function of L. A 
functional form for g(£) can then be derived by integrating 
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g(£,L) over L: 



3.2 Four selected model clusters 



Therefore, 



Lg(£)dL = g(£)L 2 max /2 



and 



(21) 



(22) 



where the expression is valid for all £ , provided L ma x is ap- 
propriately defined as the specific angular momentum per 
unit mass for a circular orbit of energy £. Given the pre- 
scription for defining r m (£), the maximal allowed angular 
momentum for a particle is easily calculated once the circu- 
lar velocity at a given radius in the potential is known. For 
the Hernquist potential one obtains, 



c(£) = rl(£) 



(24) 



With the modified definitions of g(£) and g(£,L) we 
have thus successfully adapted the equilibrium DF formal- 
ism for systems in quasi-equilibrium enabling the construc- 
tion of both f(£) and f(£,L). 



3 NUMERICAL CLUSTER MODELS 
3.1 Numerical methods 

We briefly summarize the modelling of the clusters that we 
intend to study. Initial con ditions for the cluster models are 
genera ted by means of the van de Weygaert & Bertschinger 



(1996) implementation of the Hoffman-Ribak method of con- 



strained random fields (Hoffman & Ribak 1991) in order to 
form a specific cluster at the centre of the simulation sphere. 
Each cluster model is evolve d from these constrain ed initial 
conditions by means of the Barnes & Hut 1989 treecode, 



supplemented with a galaxy formation algorithm (see van 
Kampen 1996 for details). During the evolution, galaxies are 
identified at several epochs and replaced by single 'galaxy 
particles' in order to ensure their survival, as numerical two- 
body disruption destroys alm ost all galaxies inside standard 
N-body cluster simulations (van Kampen 1995; Carlberg 



1995). Two-component models consisting of distributions of 



dark-matter background particles and galaxy particles are 
thus produced. A direct match of the simulated galaxy dis- 
tribution with the observed one sets the amplitude of the 
initial density fluctuation spectrum, and thus the present 
time in the models. The specific prescription used to define 
galaxies is necessarily somewhat arbitrary. However, we con- 
sider that it provides at least a crude representation of the 
influence galaxies have on the dark matter distribution. The 
implementation of alternate prescriptions to define 'galax- 
ies' is unlikely to substantially modify the results presented 
here. In any case, for this paper we have examined only the 
distribution and dynamics of the dominant dark component. 



We use cluster simulations from a model clu ster catalogue 
constructed by van Kampen & Katgert (1996) for the fio = 1 
standard CDM scenario. They found a consistent value 
for as (the r.m.s. density fluctuation within spheres of 
8h~ Mpc) which is used to normalise the density fluctua- 
tion spectrum. A match between models and data for the 
galaxy autocorrelation function, the richness distribution 
function, and the line-of-sight velocity dispersion fav ours 
as ~ 0.4 — .5. This is consistent with recent results of Eke 
et al. (1996 ), who find as ~ 0.5 for the same scenario. We 



therefore adopt a value of 0.5 for erg. The Hubble parameter 
Ho is set to 50 km s _1 Mpc -1 throughout this paper. 

For this paper we selected cluster Models 7, 20 and 
41 from the full catalogue. We also built the correspond- 
ing models for an fio = 0.2 CDM scenario, i.e. we kept 
the same random numbers for the initial conditions but just 
changed fio. This low-fio variant of CDM was normalized 
according to the results of Eke et al. (1996), who find from 
the number density of clusters that as scales with fio as 
0.5fio 47+0 lon °. Therefore, for fi = 0.2 we are required 
to set as ~ 1, i.e. to the unbiased value. Here we only con- 
sider the fio = 0.2 version of Model 7. These four models, 
with around 10 4 particles per cluster, form our basic simula- 
tion set. We re-simulated all models from z = 0.2 up to the 
present epoch, and saved the particle data frequently to get 
a finer timeslicing than the original run. Furthermore, since 
we do not expect to see any new galaxies forming within the 
clusters at these low redshifts, we switched off the galaxy 
formation algorithm in order to simplify the book-keeping 
of particle numbers. In Fig. 2 we show four epochs of the 
time evolution of the selected four cluster models. This fig- 
ure will be discussed in more detail in Section 5.1. 



3.3 Resolution and two-body relaxation 

The resolution of an N-body model is primarily set by the 
particle softening that has to be employed in order to re- 
duce numerical two-body relaxation effects, i.e., to make 
the simulations as collisionless as possible. However, one still 
wants to model the density distribution with sufficiently high 
spatial resolution, so a trade-off between resolution and ac- 
ceptable two-body relaxation has to be made. With this in 
mind, previous experience (van Kampen 1995) shows that 
a Plummer softening parameter of 40 kpc is a good choice 
for the dark matter particles. This corresponds to l/25th of 
the initial mean particle separation, which is well within the 
rang e of values argued for by var i ous other users of treecodes 
e.g. Bouchet fc Hernquist 198S ; Barnes & Hut 1989; Hern 



quist fc Barnes 1990} pycr fc Ip 1993| ). The softening of the 
galaxy particles depends on the half-mass radius of the orig- 
inal simulation group that was 'transformed' into a galaxy. 
Typical values fall in the range 20-50 kpc (see van Kampen 
1996), which are also well within the range usually found 
for treecodes. However, galaxy particles can have masses up 
to two hundred times larger than that of the dark matter 
particles, which means that an interaction between a galaxy 
and a dark matter particle can cause a significant change in 
the velocity of the latter, despite the softening. On the other 
hand, the dark matter particles outnumber the galaxies 50 
to 1. For the models studied in this paper the collective con- 
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tribution of galaxy particles to the relaxation of dark matter 
particles is up to five times that of the dark matter particles 
themselves (van Kampen 1996). However, for our models the 
corresponding relaxation time is still at least 6 times to, the 
present age of the universe (30io if we had just dark matter 
particles). For the chosen value of Ho here, to ~ 1 x 10 10 yr. 
So while simultaneously retaining the collisionless nature of 
the simulation, we do need to bear in mind the differing con- 
tributions of the two components to the deflection of particle 
orbits. 



4 CONSTRUCTING THE DF FOR GALAXY 
CLUSTER MODELS 

Since the positions, velocities and energies of all the particles 
that constitute a cluster model grown in a cosmological N- 
body simulation are known, the distributions h(£), g(£) and 
f(£) can be evaluated. In this section we describe the exact 
procedure followed to construct the DF for such a system. 

For any given snapshot we first compute the centre of 
mass enclosed by clumps within a 12 Mpc shell and then re- 
compute this within 2 Mpc around the most massive clump. 
Any systematic velocity with respect to the centre of mass 
is subtracted, which is equivalent to ensuring that the main 
clump defines the origin of phase-space. This in turn defines 
^t(r) and v, and thus £ through eq. (4). It is then straight- 
forward to compute the differential energy distribution h(£ ) 
as a histogram of the distribution of mass in bins of energy. 

Whereas no assumptions regarding the geometry of the 
configuration enter the definition of h(£), our aim of com- 
puting f(£) necessarily incorporates the approximation of 
spherical symmetry. This occurs during the computation of 
g(£ ) wherein although the exact potential for the particles is 
known precisely from the simulation, a spherically averaged 
value is used. The potential for Model 20 (z = 0) is shown 
in Fig. 3a. Rather than simply averaging the observed po- 
tential in each radial bin we fit an analytical form to the 
data points. Overplotted are fits to the Hernquist potential 
and the NFW potential. Both yield acceptable fits for the 
spherically averaged potential and density profile (Fig. 3b) 
within the virial radius. Note here from Fig. 3a that at the 
virial radius the potential is roughly half of its central value, 
which is still deep inside the potential well. Therefore mat- 
ter from the infall region does not significantly influence the 
dynamics within the virialized region. 

Examining several morphologically distinct cluster 
models we find that the quality of the fit to these analytic 
models depends crucially on the degree of substructure in 
the cluster which shows up as the spread in ^(r) in Fig. 3a. 
The fits get progressively worse for clumpy models since for 
these the centre of the potential is not well-determined and 
the shape departs maximally from sphericity even within 
the virial radius. 

The advantage of using the Hernquist potential fit is 
that the modified density of states, 



g(£) = 167T- 



,(£) 



(25) 
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Figure 3. The density profile and potential for Model 20 for the 
snap-shot taken at z = are plotted above. In the top panel 
(a) the points are the exact potential from the simulation and 
overplotted are the fit to the Hernquist profile (solid curve) [with 
fit parameters: * = 5.6 xlO 6 km 2 s~ 2 ; s = 2.89 Mpc] and the 
NFW profile (dashed curve) [with fit parameters V&nfw = 6.44 
XlO 6 km 2 s -2 ; s = 1.09 Mpc], The lower panel (b) shows the 
spherically averaged density profile overlaid with the Hernquist 
model (solid curve) [with fit parameters: po = 50.70 p^; s = 3.18 
Mpc] and NFW model (dashed curve) [with fit parameters: po = 
51.33 p;,; s = 1.82 Mpc], In both panels the large tick mark on 
the x-axis marks the value of riso = r max = 2.24 Mpc. 



our initial assumption of spherical sy mmetry. For these r ea- 
sons we have used the Hernquist fit (Cole & Lacey 1996) in 



can be computed analytically (see Appendix A). The dis- 
advantage of using such a fit to a given restricted form is 
in fact rather small since a much more severe limitation is 



our computations of the density of states. However, we must 
emphasize here that the entire procedure is independent of 
the precise choice of functional form for the fitted potential 
and we only use the Hernquist fit to the potential for conve- 
nience and because it does provide an elegant analytic form 
for g{£). 

The distribution function f(£) is then constructed using 
the definition in equation (9). For the snapshot at z = for 
Model 20, the computed h{£), g{£) and f(£) are shown in 
Fig. 4. From this figure we see that h(£ ) (top panel) has a 
broad distribution and has a small fraction of positive energy 
particles. The effect of our prescription for truncating the 
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Figure 4. The computed quantities for a slice of Model 20 
(Oo = 1) at z = 0. The top panel shows the differential energy 
distribution h{€) for particles inside the virial radius, plotted as a 
histogram in energy bins. The middle panel shows the density of 
states in each energy hypersurface g(£). The dashed curve is the 
usual equilibrium definition of g(£) (eq. 8) computed analytically 
using a Hernquist fit to the potential, which diverges at £ = 
(cf. Fig. 3). The solid curve uses equation (25), i.e., with the de- 
vised truncation (19) prescription presented in Fig. 1. Finally, the 
bottom panel shows the distribution function /(£) obtained by 
dividing h{£) (top panel) with the modified g{£) (solid curve in 
middle panel). The overplotted dashed line is the best particle- 
number weighted exponential fit to the curve. 



system in the computation of g(£ ) is shown in the middle 
panel, where the solid curve corresponds to the restricted 
g(£) and the dashed curve shows the asymptotic behavior 
(divergent) in the vicinity of £ = 0. In the bottom panel, 
we see that the computed /(£) is remarkably well- fit by an 
exponential function of the form f(£) = /oe _£ ^ CTfit (the 
overplotted dashed curve), in the region dominated by the 
tightly bound particles, except in the deep potential region. 
We discuss this effect further in Section 5.2. 

Defining the dimensionless central potential, 

Wo = (26) 

the equivalent r.m.s. ID velocity dispersion aiu(r) can be 
computed analytically given am and \&o by comp aring with 



the King mode l (the lowered exponential, see Binney & 



7W), 



Tremaine (1987) for details) which also has an exponential 



SW~ 5/2 aL 



with 



5 = exp(Wo) erf(VWo) 



2W \ 
3 J ' 



(27) 



(28) 



The velocity dispersion computed from the fit to the con- 
structed isotropic distribution function f(£) agrees well (to 
within 5%) with the mean ID velocity dispersion inside ri8o 
computed directly from the simulation. 

The DF can also be compared to a characteristic phase- 
space density 

3 M(r max ) 



/ 



(29) 



motivated by the expression for the central phase-space 
density of the isothermal sphere, f c = (2n)~ 3 ^ 2 p c a~ 3 . For 
Model 20, which has M(r max ) = 5.87x 10 14 M Q , r max = 2.24 
Mpc and cr m = 970 km s _1 , we find / » 1200 M Mpc~ 3 
km -3 s 3 . This fits in nicely with the distribution function 
f{£) as plotted in Fig. 4. These comparisons clearly demon- 
strate that the entire procedure produces sensible results, 
both quantitatively and qualitatively. 

To compute the anisotropic distribution function 
f{£,L) we use the definitions in equations (12), (23) and 
(24), once again with r max = ri8o, the expected virial ra- 
dius. The surface plots for h(£,L), g(£,L) and f(£,L) are 
shown in Fig. 5. The h(£,L) distribution has a coherent 
pea k at low valu es of angular momentum and energy. Un- 
like Voglis (1994 ), who analysed galaxies, we do not find a 



differentiated core and halo population for a galaxy cluster 
on equivalent scales, although there exists a small fraction 
of positive energy particles. The anisotropic DF is sharply 
peaked and falls off to zero for large £ , similar to f(£ ). To 
conclude this section, we have demonstrated that DFs can 
be constructed self-consistently for N-body model clusters 
and that sensible results are obtained. 



5 RESULTS 

In this section we analyse the dynamical states of the model 
clusters (restricted to the virialized central region by con- 
struction) using the time evolution of the DF in conjunction 
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Figure 6. The fractional change in /(£), Nf(£) during a typical 
dynamical time for the four clusters computed at z = are plotted 
above. The solid curve corresponds to Model 7 (evolved in an 
Qo = 0.2 universe), the dotted curve to Model 7, the dashed 
curve to Model 20 and the dot-dashed curve to Model 41 (all three 
evolved in an S7o = 1.0 model). The horizontal lines correspond 
to Nf = 1 and -1. 




Figure 5. The anisotropic distribution function. The figures show 
(from the top to bottom) the differential distribution h(£, L), the 
modified density of states g(£, L) and the constructed distribution 
function /(£, L) for Model 20. Here h(£, L) is in units of Mq 
Mpc -1 km -3 s 3 , g(£, L) in units of M Mpc 2 and f(£ , L) in units 
of Mq Mpc -3 km' 3 s 3 . Here £ is in units of 10 6 km 2 s~ 2 and L 
in units of Mpc km s -1 . 



with other physically interesting quantities. Before embark- 
ing on a discussion of the various distribution functions, it 
is instructive to briefly outline the evolution as seen in the 
N-body simulations for the studied models shown in Fig. 2. 

Model 7 evolving in the f2 = 0.2 universe, is relatively 
isolated and its evolution is being traced during a passively 
evolving phase. Model 7 (Qo = 1) has a significant amount 
of substructure and undergoes a violent merging event at 
the present epoch, while Model 41, which is a significantly 



more massive cluster, is being followed during one of its more 
quiet evolutionary stages. 

The three Qo — 1 models are typical for the standard 
CDM scenario in the sense that a cluster does not form 
from a single collapse but rather from a heterogenous col- 
lapse, typically from two or more diffuse sub-clumps that 
grow individually and merge to form the final cluster. Be- 
fore and after these mergers the cluster grows by accreting 
material along several filaments, which join at the position 
of the cluster. This 'secondary' infall proceeds in a somewhat 
lumpy fashion. The complex aggregation process implies 
that present-day cluster cores are often obscured by sur- 
rounding clumps. The fact that these clumps are falling to- 
wards the cluster makes it hard to distinguish them when ob- 
serving along the line of sight in redshift space. The Qo = 0.2 
cluster model is also representative for its cosmology, in the 
sense that clusters form early in this scenario, show only lit- 
tle secondary infall, and few merging clumps. We therefore 
expect this model to be relatively close to equilibrium, even 
though small clumps are still falling in and the dynamical 
timescale is about a sixth of a Hubble time. 

In the following subsections we briefly discuss what can 
be learnt from this new DF approach. The emphasis will be 
on demonstrating the usefulness and reliablitity of our con- 
struction of distribution functions, using the N-body models 
to illustrate the results. Detailed discussions of the dynamics 
of cluster formation and the cosmological implications using 
this approach are deferred to later papers. 



5.1 Time scales and quasi-equilibrium 

As anticipated in Sections 2.1 and 2.3, it is evident from the 
evolution plots in Figs. 7, 8 and 9, that in general, even over 
relatively short time scales (compared to the Hubble time), 
the dynamical state changes perceptibly, thus validating the 
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quasi-equilibrium approach of treating the system as one 
evolving over two characteristic but disparate time-scales. 

To test the assumption that clusters are indeed in a 
quasi-equilibrium state at least in the interim between major 
merger events, we calculate the fractional change in f(£) 
during one dyn amical time idyn = y 37r/16Gp (Binney & 
Tremai le 1987), which is equal to 0.25io (for the flat model), 
as p = piso = 180pb for the virialized region under study (by 
definition) . At z — we compute the dimensionless quantity 
Nf(£) defined to be a measure of the stationarity, 



Nf(£) = idyn( 



dt 



(30) 



plotted in Fig. 6. In order to get a quantitative measure 
for the stationarity of /(£), we calculate the h(£ )-weighted 
averages of Nf(£) for all four models and hence define the 
'stationarity index' n/ as, 



/ \N f (£)\h(£)d£ 



M 



(31) 



where M is the total mass enclosed within riso for the sys- 
tem. The values of n/ that we find are: 0.32, 3.3, 3.7, and 0.96 
for Model 7 (fi = 0.2), 7 (fl = 1), 20, and 41 respectively. 
We interpret this as follows: only Model 7 for the Qo = 0.2 
case can really be considered to be in quasi-equilibrium, with 
clear skewing evolution of /(f), while the same model within 
an fio = 1 universe, where it is undergoing a major merger, 
and Model 20 are furthest from an equilibrium state. Model 
41 is a 'critical' case: there certainly is on-going secondary 
infall, but the cluster is almost massive enough for the viri- 
alized region to remain in quasi-equilibrium. Therefore, we 
roughly delineate three distinct regimes on the basis of the 
values of the 'stationarity' index nf. nf < 1 corresponds 
to the regime where the quasi-equilibrium description is an 
adequate one, nf = 1 to the critical case and nf > 1 to the 
case when the system has departed from quasi-equilibrium 
while recuperating from a merger. 



5.2 Signatures of substructure 

As mentioned in Section 4 there appears to be a central 
peak in f(£) for most of the models. The peak is partic- 
ularly marked when there is significant substructure inside 
the studied region. We therefore associate this peak with 
infalling sub-clumps outside the cluster core which already 
have a deep 'local' central potential, but contain an insignif- 
icant fraction of the mass of the system. Indeed, the effect is 
not seen in h(£ ) confirming the above picture. The central 
peaks in f(£) are thus partially an artifact of the assump- 
tion of spherical symmetry of the potential, but are also a 
signature of transient substructure. Once the clumps are in- 
corporated or fall in to the centre of the cluster, the peak 
disappears. 

The radial profile for the velocity anisotropy parameter, 



0(r) 



1 



<r r 2 (r)' 



(32) 



is correlated to the presence of substructure (see Fig. 2). In 
the absence of significant sub-clustering the velocity distri- 
bution is found to be roughly isotropic but erratic (f3 « 0) 
and for quiescent evolution j3 > 0, indicating the presence 



of primarily radial orbits, which are most likely to originate 
in the infall region. 

The presence of substructure in the cluster models being 
progressively subsumed into the main clump can be seen 
in the (significant) evolution of the ID velocity dispersion 
profiles as well. For example, consequent to a merger the 
mean value increases as shown in the top left panel of Fig. 8. 
Furthermore, the overall shape seems to be well correlated 
with substructure and secondary infall. This promises to be 
one of the important applications of the DF formalism. 



5.3 Evolution in the density profiles 

Comparatively little systematic evolution is seen in the den- 
sity profiles, aside from the progressive worsening of the fits 
to the analytic form occuring due to the presence of sub- 
clumps causing significant departures from sphericity, and 
small changes in the value of the central density. Density 
profiles seem to be set up rather early in the evolutionary 
history of the cluster. 



5.4 The central potential and degree of central 
condensation 

The dimensionless depth of the potential, Wo as defined 
in equation (26), is a measure of the central concentration 
of the system. For the models studied, we find Wo ~ 7 
for the flat models and Wo w 10.5 for the open model. 
This is comparable to what is found for elliptical galax- 
ies, where 8 <> W & 11 ( |Hjorth fc Madsen 1995] ), and 
uncollapsed globular clusters for which Wo 5s 10 (eg. Bin- 
ney & Tremaine 1987). Clusters evolving in an Einstein-de 
Sitter universe have comparatively shallower central poten- 
tials since the system is constantly disturbed by infalling 
sub-clumps and is hence unable to build up a concentrated 
core in the elapsed time. By contrast, a cluster evolving in 
an open universe is comparatively more homogeneous and 
relatively undisturbed. It evolves quiescently essentially due 
to the fact that dark matter particles are allowed to plunge 
down into the centre on radial orbits (as evidenced by j3 > 
seen in the middle left panel of Fig. 7) hence building up to 
high central densities. 

The observed trends in the central concentration can 
also be potentially related to features of the initial condi- 
tions by estimating the characteristic phase-space density 
/. The typical maximum phase-space density is about 100 
times larger than this characteristic value, similar to what 
has been found for elliptical galaxies. Since the maximum 
phase-space density cannot in crease during collisionless evo- 



lution for an isolated system (Tremaine, Henon, & Lynden 



Bell 1986 ), the computed value of / allows us to under- 
stand the relative efficiency required of any relaxation pro- 
cess, given characteristic time-scales of mediation, in order 
to produce the present configuration from the inhomoge- 
neous initial cosmological conditions from which a cluster 
first assembles. 



5.5 Energy distribution 

There is a clear difference between the differential energy dis- 
tributions, h(£), of the open and flat models. The energy dis- 
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Figure 7. The evolution of computed quantities for a slice of Model 7 (do = 0.2) for four redshifts chosen to correspond to constant 
time intervals. As a function of physical radius, the left panel shows the measured one-dimensional velocity dispersion <r 1D (r) (top), the 
anisotropy parameter /3(r) (as defined in eq. 32) (middle), and the density p(r) relative to the background density p b (bottom). In the 
right panel the differential energy distribution h(£) (top) and the computed isotropic distribution function f(£) (middle) are plotted as 
a function of the binding energy per unit mass £. 
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Figure 8. Same as Fig. 7 for Model 7 within an S7q = 1 Universe. 



tribution for the open model (Fig. 7) is very similar to that 
of elliptical galaxies where a 'negative temperature' expo- 
nential function, cut-off close to the escape energy {£ = 0), 
provides a good fit the h(£ ) (Binney 1982), i.e., most of the 
particles are loosely bound. By contrast, the energy distri- 
bution of the flat models are better approximated by gaus- 



sians with a tail of positive-energy particles. The advantage 
of our adapted DF formalism is that it allows the inclusion 
of positive energy particles as well. For the galaxy cluster 
models evolved in a flat universe we find that the particles 
with positive energies that are located preferentially in the 
centre of the cluster, despite the deep potential well. Energy 



14 Priyamvada Natarajan, Jens Hjorth and Eelco van Kampen 



1500 - 



E 1000 




500 




1.0 



z= 0.00: o obs =1 374, a fit =1 1 63 




-0.5 



-1.0 L 
0.0 



10" 



10 c 



* 10 7 



10 c 



10 



10° 
10 



10 c 



10 L 



'& 10 4 



- 10' 



10 L 



10"' 





10 6 




10 5 




10 4 


CL 


10 3 


i— 




Q_ 


10 2 




10 1 




10° 




10" 1 



0.5 1.0 1.5 2.0 2.5 3.0 
r [Mpc] 



z= 0.00: p =179.1p b , s=1.78 




0.1 



1.0 
r [Mpc] 



10.0 





z= 0.00: W = 7.19, a, =11 82 



5 10 

e [10 6 km 2 s 2 ] 



z= 0.20 
z= 0.13 
z= 0.06 
z= 0.00 



Model 41, Q =1.0 



Figure 9. Same as Fig. 7 for Model 41 (fto = !)• 



exchange due to violent merger events as seen in Model 7 
for flo = 1-0, for instance, are a potential source for gener- 
ating these unbound particles, as we clearly see a significant 
fraction of them just before and during a merger. We also 
detect simultaneously a surge in the maximum value of £ , 
indicating significant energy diffusion. We hope to address 



this issue in the context of examining the relaxation process 
in more detail in a future paper. 
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6 SUMMARY AND CONCLUSIONS 

In this paper we have applied an adapted version of the 
distribution function formalism of equilibrium stellar dy- 
namics to numerical models of clusters of galaxies. Clusters 
are generally not in equilibrium as on-going secondary infall 
and the relatively long dynamical time-scale prevent them 
from reaching an equilibrium configuration. However, as ev- 
idenced from observations of clusters, this deviation from 
equilibrium is not expected to be too severe, at least within 
the finite central virialized region, and so a quasi-equilibrium 
description can be formulated enabling the construction of 
the distribution function for the N-body cluster models. 

Isotropic /(£) and anisotropic f(£,L) distribution func- 
tions can be computed using the differential energy distri- 
bution and the density of available states in phase space for 
the system. The accessible phase space however, needs to be 
restricted in a physically meaningful and self-consistent fash- 
ion. We invoke the finite elapsed time since formation and 
space available as providing a natural way to limit phase 
space. We argue that this prescription can be suitably ap- 
plied when the analysis is restricted to the virialized region, 
which is defined by r < rigo. This is motivated by the fact 
that the particles that once enter and get bound within this 
region are not likely to leave. Even particles that do be- 
come unbound due to energy-exchange processes are likely 
to lose energy fairly rapidly in the dense cluster environment 
on very short time-scales, i.e., their occupancy of states in 
phase space is automatically limited by the fact that parti- 
cles are expected to change energy again within a relatively 
short time. It is precisely the existence of these two dis- 
parate time scales (the long dynamical time scale and the 
short characteristic time scale for energy exchange) which 
implies that the evolution of the system can be effectively 
decoupled into 'fast' and 'slow' components and hence al- 
lows the description of quasi-equilibrium states somewhat 
within the equilibrium formalism. In the anisotropic case it 
is additionally assumed that g(£,L) is separable. 

We show that this construction can be successfully ac- 
complished for N-body cluster models taken from a large 
sample of models for the standard CDM cosmological sce- 
nario, and that the constructed DFs do show similarities to 
equilibrium ones. The velocity dispersion as obtained from 
an exponential fit to f(£) agrees well with that obtained 
directly from the particle distribution in the N-body simu- 
lation (agrees to within 5 % typically). This is not only true 
for specific redshift slices, but also during the course of evo- 
lution. The only exceptions are during strong merger events, 
whence the difference can be a few hundred km s -1 . The 
characteristic phase-space density, calculated from global 
properties of the particle distribution in the simulation and 
from the computed DF also agree well. 

Furthermore, we show by following the time variation 
of the DFs that they are good tracers of dynamical events 
like cluster mergers and infall of clumps. There is notable 
evolution with redshift for most models, indicating on-going 
dynamical evolution. The exception is the cluster model sim- 
ulated in an fio = 0.2 CDM universe, which shows only 
mild evolution and can be said to be in quasi-equilibrium. 
It is also the only model to show a linear rise of the velocity 
anisotropy parameter j3 from zero in the centre of the cluster 
to about 0.5 at riso- This cluster bears the closest resem- 



blance to an elliptical galaxy, further demonstrated by the 
shape of h(£ ) and the degree of central concentration which 
is well within the range typically estimated for ellipticals. 
We also find that the density profiles within the virialized 
region show practically no evolution and seem to be set up 
quite early in the dynamical history. 

As to the possibilities of using the constructed DF to de- 
scribe the dynamical state of a galaxy cluster, it can clearly 
distinguish the merging clusters from the quiet ones, as well 
as the ones dominated by secondary infall. The evolution 
of the DF, even over a relatively short time, reveals the 
on-going kinematic activity. The comparison of global prop- 
erties like the velocity dispersion as obtained from a fit to 
/(£) to those calculated directly from the simulation pro- 
vide a new indicator to establish whether the cluster is close 
to equilibrium. 

In a recent study of the evolution of cluster-scale dark 
halos, Formen, Bouchet, & White (1996) also characterize 



the halo formation process as being composed of alternat- 
ing merging and relaxation phases. The general picture that 
emerges from our analysis is one of clusters evolving qui- 
escently with on-going collective energy exchange (which 
we define as quasi-equilibrium states) punctuated by violent 
merging phases followed by fairly brief, but yet well-defined 
recovery phases. During the passively evolving phases, we 
find that the DF f(£) is well-approximated by an exponen- 
tial and significant deviations occur only during the merging 
and recouping phases. In this treatment, we have not quan- 
tified the frequency or nature of energy exchange processes, 
aside from the observation that they do indeed occur. Be- 
sides possibly allowing a better understanding of the physics 
of energy exchange, this approach also potentially affords 
new discriminants of the underlying cosmological model. 

In conclusion, we have shown that distribution functions 
can be successfully constructed for the dark matter com- 
ponent of evolving, non-isolated clusters of galaxies grown 
in N-body simulations. The necessary machinery is thus in 
place for future work in which we hope to address the use 
of distribution functions as a new probe of cluster dynamics 
and evolution, fundamental cosmological parameters, and 
the possible origin of 'universal' cluster density profiles. 
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APPENDIX A: TRUNCATED DENSITIES OF STATE FOR VARIOUS POTENTIALS 



For the Hernquist (1990) potential used in this paper, 
fro 



*(r) = 



1 + 



where s is a characteristic scale radius, the truncated density of states can be expressed analytically as, 

g{£) Vfrp - £ (8£ 2 + 10£ fro - 3 fr 2 ) s 3 
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Alternate popular potentials do not completely reduce to analytic expressions. For the NFW potential (Navarro, Frenk 
& White 19961), 



>Tf , s *NFWS , . r> 

fr(r) = ln(l + -), 

r s 

the truncated density of states is 
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For the Jaffe model ( |jaffe 1983| ), 

*(r) = frj ln(l + — ), 
r 

the truncated density of states is 



g(£) _ *jO 



„ , - , ■ ^ + m % /fr J ln(l + ^)-£. 
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For the singular isothermal sphere (Binney & Tremaine 1987), 

*(r) = */ ln(i), 

the truncated density of states is 
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